The genetic basis of conspicuous coloration in the Guadeloupean anole: Evolution by sexual and ecological selection

Abstract Understanding how natural selection acts on the genome and contributes to the process of speciation is a primary aim of the study of evolution. Here we used natural variation in two subspecies of the Guadeloupean anole (Anolis marmoratus ssp.), from the island of Guadeloupe in the Lesser Antilles, to explore the genomic basis of adaptation and speciation in Anolis lizards. These subspecies inhabit distinct ecological environments and display marked differences in adult male color and pattern. We sequenced the complete genomes of 20 anoles, 10 from each subspecies, at 1.4× coverage. We used genome‐wide scans of population differentiation, allele frequency spectrum, and linkage disequilibrium to characterize the genomic architecture within and between the subspecies. While most of the genome was undifferentiated, we observed five large divergent regions. Within these regions we identified blocks, 5 kb pairs in length, enriched for fixed single nucleotide polymorphisms. These blocks encompass 97 genes, two of which are candidate pigmentation genes. One is melanophilin (mlph), which helps transport melanosomes within melanocytes. The other is a cluster of differentiation 36 (cd36), which regulates carotenoid pigment sequestration. We used high‐pressure liquid chromatography to confirm that carotenoid pigments are significantly more abundant in the conspicuous orange‐pigmented skin of male A. m. marmoratus suggesting that cd36 may be regulating pigment deposition in this tissue. We identified for the first time a carotenoid gene that is a potential target of divergent sexual selection and may be contributing to the early stages of speciation in Anolis lizards.


| INTRODUC TI ON
Characterizing the genetic architecture of speciation is fundamental to understanding how new species evolve. While the mechanisms that drive speciation, such as natural selection, sexual selection, and genetic drift are well understood, there are few examples of how these processes act on the genome (Ellegren et al., 2012;Jones et al., 2012;Kronforst et al., 2013;Poelstra et al., 2014). Speciation research has traditionally focused on the genes that underlie reproductive isolation and not on those genes that contribute to the earliest stages of species divergence (Coyne & Orr, 2004;Endler, 1990;Moyle & Payseur, 2009;Orr et al., 2004;Presgraves, 2010). It can be argued that genes that initiate divergence are of greater interest to evolutionary biologists because reproductive isolation typically arises as a by-product of selection, drift, or founder effects in isolated populations (Andersson, 1994;Kolbe et al., 2012;Ritchie, 2007;Shaw & Mullen, 2011). Thus, characterizing the genetic architecture of divergence allows for the identification of the genes and genomic loci vital to initiate the process of speciation.
Speciation can occur with varying degrees of connectivity between populations (Schluter, 2001) with new species arising in a continuum from sympatry, when the ranges of populations overlap entirely, to allopatry when populations are completely separated (Butlin et al., 2008). Here we focus on parapatric speciation, which occurs when populations diverge, but continue to overlap and interbreed for part of their range (Endler, 1977). Parapatric speciation is ubiquitous and is often driven by diversifying ecological selection along an environmental gradient (Doebeli & Dieckmann, 2003;Fu & Li, 1993). While divergent ecological selection is thought to play a central role in parapatric speciation (Kondrashov & Kondrashov, 1999;Schluter, 2009;Schneider, 1996), sexual selection also may contribute to divergence if the sexually selected trait is differentially detectable along the environmental gradient (Lande, 1982;Maan & Seehausen, 2011). Thus, populations in the earliest stages of parapatric speciation are expected to show divergence at genes and loci that contribute to phenotypes under ecological selection, sexual selection, or both processes.
Our goal was to characterize the genomic architecture of ecological diversification in two subspecies of Anolis lizards. Anoles are a model organism for the study of speciation, and represent one of the largest adaptive radiations of any terrestrial vertebrate (Losos, 2009).
In the Caribbean, anoles have evolved similar body plans when species independently adapted to almost identical environments on different islands in what is considered a classic example of convergent evolution (Turner et al., 2005;Williams, 1972). Anoles also display a remarkable variety of colors and patterns both as adaptations for crypsis and for species and mate recognition. Most species of anoles possess a brightly colored extensible throat fan called a dewlap, which is used in inter-and intraspecific communication (Nicholson et al., 2007;Ord & Martins, 2006;Vanhooydonck et al., 2009). The differences in dewlap color and or pattern plays an important role restricting gene flow between species that is habitat dependent (Ng et al., 2017). However, in the early stages of speciation, anoles typically diverge first in male dorsal color and pattern (Losos, 2009), but see Stapley et al. (2011). This suggests that variation in dorsal coloration might arise as a result of natural selection, as well as female mate choice, although this has yet to be demonstrated. Dorsal coloration and pattern is shaped by both ecological and sexual selection as a consequence of local adaptation to spatially varying environments (Medina et al., 2017;Muñoz et al., 2013). Thus, divergence in dorsal coloration may initiate ecological divergence and be an important driver in speciation and adaptive radiation in anoles.
The bright colors of lizards are known to consist of two types of pigments: carotenoids and pteridines (reviewed in Olsson et al., 2013). Both carotenoid and pteridine pigments have been identified in the red and orange dewlaps of anoles (Macedonia et al., 2000). Although they produce similar colors, carotenoid pigments are sequestered from dietary sources (Olson & Owens, 1998) while pteridine pigments are synthesized endogenously from purines (i.e., guanosine-5′-triphosphate; Ziegler et al., 2000). Carotenoid pigments are of particular interest because they often color the bright integument of vertebrates (Goodwin, 1984), play important roles in immune function (Saks et al., 2003), and have antioxidant properties (Stahl & Sies, 2003). Thus, carotenoids are potentially honest indicators of quality if they color a sexually selected trait (Endler, 1980(Endler, , 1983Kodric-Brown, 1985;Olson & Owens, 1998). However, carotenoid research on anole dewlap color has found no influence of food access or carotenoid supplementation (Steffen et al., 2010), but the role of carotenoids on anole dorsal color is unknown.
Although 11 candidate genes have been identified as potentially associated with carotenoid pigmentation (Walsh et al., 2012), the pathway that regulates carotenoid sequestration and deposition in vertebrates remains largely uncharacterized. Conversely, the genetic basis of pteridine pigment synthesis has been characterized in zebrafish and fruit flies (Dupont, 1958;Forrest & Mitchell, 1954), and is regulated by at least seven genes (Braasch et al., 2007;Kim et al., 2013;Ziegler et al., 2000). In anoles, De Mello et al. (2021) recently identified 13 candidate genes involved in anole dewlap color and color pattern. However, none of the genes in either pathway have been investigated in anole dorsal coloration. Because conspicuous coloration is one of the first traits to diverge in anoles, it is likely that either carotenoid or pteridine genes play a role in the speciation of anoles.
On the island of Guadeloupe in the Lesser Antilles, the Leopard or Guadeloupean anole (Anolis marmoratus) represents a unique example of an anole species in the early stages of parapatric speciation (Lazell, 1963(Lazell, , 1972Muñoz et al., 2013;Schneider, 1996).
Anolis mamoratus subspecies are primarily defined by differences in male color and pattern and are found in a variety of habitats ranging from wet cloud forest to arid scrub forest (Lazell, 1963). The two main islands are inhabited by six subspecies: five on Basse-Terre and two on Grande-Terre with one subspecies found on both islands (Lazell, 1963). Pairwise comparisons of mitochondrial DNA reveal very limited geographic structure between the subspecies (Schneider, 1996), and microsatellite loci genotyped in 11 populations on Grande-Terre reveal no evidence of isolation by distance (Muñoz et al., 2013). A population genomic analysis also found high rates of gene flow among populations on both main islands (T. J. McGreevy, N. G. Crawford, P. Legreneur, & C. J. Schneider, unpublished data). This lack of genetic differentiation, combined with substantial phenotypic differences in male coloration, provides a unique opportunity to characterize the genetic architecture of phenotypic divergence. Importantly, because visual signals are key components of species recognition in anoles (Losos, 2009;Macedonia et al., 2013), our analysis is relevant to understanding the early stages of speciation in anoles.
Here we focus on two subspecies of A. marmoratus: A. m. speciosus and A. m unpublished data). Anolis m. speciosus is found in mesic woodlands, and males have a conspicuous blue wash on their heads that is conspicuous in the light environment characteristic of their habitat (Muñoz et al., 2013). The blue coloration is consistent with expectations for effective signaling colors in large gaps in forests which are enriched for blue wavelengths of light (Endler, 1993). In contrast, Anolis m. marmoratus inhabits the southeastern coast of Basse-Terre, and is found in closed tropical rainforest inhabiting what would be considered small gap light environments (sensu Endler, 1993).
Males have conspicuous orange marbling on their heads consistent with expectations for conspicuous visual signals in small gap light environments which are enriched in the red-orange wavelengths (Endler, 1993). There is clinal variation in head color across the range of the two subspecies with populations in the central eastern lowlands of Basse-Terre displaying mixed phenotypes. The amount of orange on the head increases southward from just north of the town of Goyave and reaches its peak in the vicinity of Capesterre.
Inversely, the amount of blue increases (and orange decreases) as one moves northward from south of Goyave to northeastern Basse-Terre and southeastern Grande-Terre, where blue coloration peaks in the vicinity of Pointe-a-Pitre. The conspicuously different phenotypes that we sampled in this study represent the end points of a cline in head coloration which exists despite the levels of gene flow sufficient to homogenize most of the genome. Our working hypothesis is that the conspicuous differences in head coloration result from divergent sexual selection (both intra-and inter-sexual selection) for effective signals driven by differences in light environment.
Using measures of light environment, skin reflectance, highperformance liquid chromatography (HPLC) of carotenoid pigments, double digest restriction-site-associated DNA sequencing (ddRADseq), and low-coverage whole genome sequencing we address the following questions in regard to the early stages of speciation in anoles: (1) which phenotypic traits are under divergent selection; (2) what is the genomic basis of divergence in those traits; and, because divergence involves orange coloration, (3) are genes associated with pteridine or carotenoid pigmentation involved? We hypothesize that due to the high rates of gene flow between A. m. speciosus and A. m. marmoratus, there will be small islands of divergence that are associated with genes involved in phenotypic divergence. And, because F I G U R E 1 Representative photos of adult male Anolis marmoratus: (a) Anolis m. marmoratus from Capesterre and (b) Anolis m. speciosus from Pointe-à-Pitre. (c) depicts the approximate ranges of each subspecies with the color gradient representing the approximate transition between their dorsolateral head colors, the two collecting localities, and the clinal variation between them. The insert displays the location of the island of Guadeloupe relative to other islands in the Lesser Antilles.
coloration plays an important role in speciation in the adaptive radiation of Caribbean Anolis, our findings will be relevant for understanding color evolution across the radiation. A white barium sulfate standard was used to zero the instrument.

| Sample collection and phenotyping
Special care was taken to ensure spectral measurements were not taken when the animals were stressed and in a darkened phenotype color. Reflectance measures were converted to hue following Endler (1990).  (Table S1). Each sample had their genomic DNA standardized to 500 ng and was simultaneously digested with MluCl and NlaIII at 4°C overnight. Each sample was uniquely barcoded with adapters that ranged in size from 34 to 53 nucleotides. Samples were pooled at equal concentrations and subsequently size selected between 182 and 222 bp with a Pippin

| Whole genome library preparation and sequencing
Prep® electrophoresis system (Sage Science). Samples were singleread sequenced on one lane of 50 bp Illumina HiSeq™ 2000 sequencing system.

| Quality control, read alignment, and SNP/ indel calling
Reads were aligned to the A. carolinensis reference genome (AnoCar2, ensembl version 67) with Stampy (1.0.18, r1526). Stampy can incorporate prior information about expected sequence divergence from the reference genome, so we set the "substitution rate" option to 0.13 and ran Stampy with the "--sensitive" flag set.
The Broad Institute's Genome Analysis Toolkit (GATK version 2.2-16, g9f648cb) was used to identify single nucleotide variants (SNVs) and short insertions and deletions (indels). Following GATK's best practices, duplicate reads were identified with "MarkDuplicates" with Picard command line tools (version 1.72). Then indels were identified using GATK's "RealignerTargetCreator" and reads were locally realigned using "IndelRealigner". SNVs and short indels were called using "UnifiedGenotyper" with permissive setting ("minimum base quality scores" were set to 2). Variants were then recalibrated using high-quality SNVs identified in ddRAD dataset (SNV quality >500). Because we were concerned that the ddRAD dataset did not contain enough indels to provide an accurate truth set (~ 3 indels/ thousand SNVs), we simply discarded those indels whose quality score fell into the lowest quantile scores.
Because BreakDancerMax identification requires high coverage to accurately identify structural variant, we applied the algorithm only at the level of populations (~14× coverage). We also only analyzed reads with mapping quality scores of at least five (i.e., reasonably likely to be uniquely mapped).
To calculate G ′′ ST (Meirmans & Hedrick, 2011) and to identify fixed SNPs, variant call format (VCF) files were processed with custom python code (https://github.com/ngcra wford/ pypgen). We used ANNOVAR (Wang et al., 2010) to annotate SNVs that intersect with genes as well as to identify synonymous and non-synonymous SNVs. Weir and Cockerham's F ST (Weir & Cockerham, 1984) and Tajima's D (Tajima, 1989) for each subspecies was calculated with vcftools (Danecek et al., 2011). G ′′ ST and Tajima's D were calculated in 5 kilobase pair (kbp) nonoverlapping blocks, and only SNVs for which there were five samples per population were used in these calculations. We choose 5 kbp as our block size because this was the minimum block size where each block contained, on average, enough SNVs (x = 109.12 ± 52.00 stdev) to calculate summary statistics without sacrificing precision.
To measure the extent of linkage disequilibrium along the genome, we used Beagle (version 3.3.2) to phase our genotypes (Browning & Browning, 2007). We then used VCF tools to calculate correlation coefficients between all SNVs in 25 kbp nonoverlapping blocks. Blocks smaller than 25 kbp did not contain enough variants to accurately fit decay curves. Because r 2 is sensitive to rare alleles (Remington et al., 2001), we only included SNVs where the minor allele frequency was ≥0.2 and where the proportion of missing data was less than 20%. To each window we fitted a decay curve (Hill & Weir, 1988;Weir & Hill, 1986) and measured the fitted r 2 at the midpoint of the window (Alhaddad et al., 2013).

| Carotenoid HPLC
Skin tissue was excised from the most pigmented portion of the dorsolateral region above and behind the eye. Tissues were immediately frozen in liquid nitrogen and transferred to a −80°C freezer for longterm storage. All tissues were collected at the same time and stored in the dark.

All HPLC steps were performed at the Carotenoids and Health
Laboratory at the Tufts School of Nutrition. Samples were extracted in a darkened room under red light to reduce photodegradation of pigments. After thawing, tissues were weighed on a Mettler Toledo AX26 Comparator high sensitivity balance. Each tissue was individually processed for HPLC analysis.
To this volume, 1 mL of 0.85% saline and 100 μL of internal standard (echinenone diluted to absorb at 0.02) were added. The sample was then vortexed for 30 s.
To separate the CHCl 3 :MeOH mixture, the sample was centrifuged for 10 min at 800 g at 4°C. The lower CHCl 3 layer was transferred to a fresh test tube and evaporated to dryness in a 40°C water bath. Additional pigments dissolved in the aqueous phase were extracted by adding 3 mL of hexane to the remaining aqueous layer, vortexed for 1 min, and centrifuged for 10 min at 800 g. The upper hexane layer was transferred to the dried sample and evaporated to dryness. The dried pigments were resuspended in 100 μL of 100% ethanol, vortexed for 1 min, sonicated for 30 s, and transferred to an autosampler vial for HPLC analysis.

| HPLC analysis
Each sample was processed on a Waters HPLC instrument with a Waters 2996 Photodiode Array Detector, a Waters 2695 Separations Module, and a Waters 474 Scanning Fluorescence detector. A YMC(tm) Carotenoid column S-3, 3.0 × 150 mm (# CT99S031503WT) was used to separate the pigments.
The running method used two solvents, A and B. Solvent A was composed of: 0.85% methanol, 0.12% methyl tert-butyl ether (MBTE), 0.03% H 2 0, and 0.45 g/L ammonium acetate. Solvent B was composed of 0.08% methanol, 0.90% MBTE, 0.02% H 2 0, and 0.2 g/L ammonium acetate. Samples were run at a constant flow of 0.4 mL/min with A and B buffers. After the column was equilibrated, 20 μL of sample was injected. Then the following separation protocol was run: 21 min linear gradient from 100% solvent A to 45% solvent A, 1 min at 45% solvent A, 11 min linear gradient to 5% solvent A, 4 min at 5% solvent A, 2 min linear gradient to 100% solvent A, and 21 min at 100% solvent A. Peaks were identified using Empower 3 Software (Waters) and external calibration standards from Sigma-Aldrich.

| Phenotypic measurements
To characterize differences in color and pattern between male A. m. marmoratus and A. m. speciosus, we measured the hue of the dewlap, dorsolateral head, skin around the eye (eye-ring), dorsolateral body, and the tail. Hue relates measures of reflectance of the discrete colors in the visible spectrum (e.g., red, yellow, blue, and green; Endler, 1990). The only significant difference in pigmentation was in dorsolateral head coloration.
Color reflectance from ~32 millimeter skin patches on the dewlap, dorsolateral head, eye-ring, dorsolateral body, and tail was measured from 10 individuals from each population (Figure 2). Reflectance was converted to hue following Endler and Mielke (2005 ; Table S2) Table S3). A Tukey test identified 23 (p < .05) significant differences, but most were between different types of skin patches and did not explain differences between subspecies. Only the mean hue between dorsolateral head patches in A. m. marmoratus and A. m. speciosus was significantly different (p < .001).
Because orange coloration may result from carotenoid pigments, we used HPLC to measure carotenoids in the head skin. We measured carotenoid pigments from five male A. m. marmoratus and five male A. m. speciosus (Figure 3; Table S4). Carotenoids were classified into five types: β-cryptoxanthin, lutein, zeaxanthin, carotenes, and esterified xanthophylls ( Figure S1 Table S5]. A Tukey test identified nine significant differences (p < .05), but most were between different types of carotenoids and not biologically relevant. Only the mean concentration of esterified xanthophyll pigments between A. m. marmoratus and A. m. speciosus was significant (p < .001).

| Sequencing and genomic variation
We sequenced 20 samples, 10 from A. m. marmoratus and 10 from A. m. speciosus, in two lanes of 101 bp paired-end Illumina reads (

| Summary statistics
We used G ′′ ST , an F-statistic similar to F ST , to measure the degree of genetic differentiation between subspecies (Table 1). Unlike F ST ,

G ′′
ST accounts for multi-allelic sites, small sample sizes, and a small number of sampled populations (Meirmans & Hedrick, 2011).

We measured G ′′
ST across the genome in 5 kbp blocks (334,498 in total). Globally, mean G ′′ ST does not vary among chromosomes (x = 0.166 ± 0.106 stdev) with the exception of microchromosome LGb (x = 0.257 ± 0.148 stdev), which had a value 1.5 times higher. We identified both the top 1% (N = 3345, x = 0.5399) and 5% of nonoverlapping blocks (N = 16,725, x = 0.3569). We defined outliers as blocks falling within the top 1% of G ′′ ST outliers. We also calculated Tajima's D, a statistic which summarizes the allele frequency spectrum (Tajima, 1989). Tajima We identified five regions of increased divergence, when the top 1% of 5 kbp G ′′ ST blocks were physically clustered together on macrochromosomes 1, 2, 3, 5, and 6 ( Figure 4). Within the clusters of outliers, Tajima's D was significantly reduced and LD was significantly increased relative to the rest of the genome at different alleles in
When we examined the distribution of SNPs fixed between populations, we observed that they generally clustered at the center of the divergent regions. These regions intersect 97 genes of which only two, cluster of differentiation 36 (cd36) and melanophilin (mlph), are associated with pigmentation (Table S8). Eight SNVs fall within mlph introns and five within cd36 introns.

| DISCUSS ION
Here we characterized dermal color differences in the adult males of two subspecies of A. marmoratus. We identified significant differences in hue and carotenoid content in the dorsolateral skin of the head. We show that genetic divergence occurs at only a few genomic regions comprising approximately 2% of their genome. In these regions, we identified two pigmentation genes: mlph, which is involved in melanosome transport, and cd36, which regulates carotenoid sequestration.
These results suggest divergent selection on coloration, in this case likely sexual selection for effective visual signals, may initiate adaptive divergence in A. marmoratus populations and may be important in color divergence associated with speciation in anoles more broadly.
However, direct tests of ecological and sexual selection are needed to determine the extent of their influence on divergence.

| Phenotypic divergence
Because the head coloration is sex-specific, conspicuous, and does not itself contribute to niche specialization, it is likely evolving primarily by sexual selection (Andersson, 1994;Ritchie, 2007). We hypothesize that local adaptation to light environment is a key to understanding divergence in color phenotypes. In particular, Endler's framework developed in his seminal 1993 paper, "The color of light in the forest and its implications" is relevant here. Endler (1993) shows that different light environments exist in the forest and suggests that visual signaling colors could be locally adapted to those light environments. We suspect that is what is happening in the anoles of Guadeloupe, which also has been seen with other anole species (e.g.,  (Muñoz et al., 2013). This pattern of general conspicuousness also was noted by Williams and Rand (1977) across species of Anolis more broadly. We hypothesize that the striking variation between the described subspecies A. m. marmoratus and A. m. speciosus represents local adaptation to different light environments at either end of a habitat and light environment cline between mesic forest and closed rainforest. Overall, the light environment, where the nominate subspecies marmoratus is found, is characteristic of light in the rainforest where it is dim overall, but small gaps generate a light environment enriched for red-orange wavelengths (Endler, 1993 Note: r 2 was calculated in 25 kb blocks. F ST was calculated following Weir and Cockerham (1984) using VCFtools (Danecek et al., 2011). G ′′ ST was calculated with custom python code. Tajima's D was calculated with VCFtools. r 2 was calculated from a combination of R and custom python code. Genomic background includes all blocks with positive F ST values excluding the top 5% of outliers.
enriched for blue wavelengths (Endler, 1993). Therefore, we suspect that local adaptation of signaling colors to light environment is driving the divergence of these populations and that only a small portion of the genome is diverging-the remainder is homogenized by the high levels of overall gene flow among the populations. Because visual cues are a critical part of the Anolis species recognition system, divergence in color phenotypes and the associated genes has direct implications for understanding the early stages in the process of speciation in anoles more broadly.

| Genomic architecture and identification of divergent regions
In both subspecies, genome-wide genetic divergence, as measured Within the top 1% of outlier loci, we observed that Tajima's D was significantly lower than the genome-wide average. Reduced Tajima's D is indicative of purifying selection which reduces genetic variation resulting in an excess of low frequency polymorphisms.
We also observed increased linkage disequilibrium in these outlier loci. Increased linkage disequilibrium is expected in regions under purifying selection because decreased recombination increases the chance that selected alleles will remain in linkage disequilibrium.
Together the regions of localized divergence combined with the reduced Tajima's D and increased linkage disequilibrium strongly suggest that selection is acting on the top 1% of divergent genomic regions.
We observed that these outliers appear to cluster into five divergent regions. Divergent regions are a common feature of speciation with gene flow and have been observed in butterflies, fish, and birds (Ellegren et al., 2012;Jones et al., 2012;Kronforst et al., 2013;Poelstra et al., 2014). In the literature these clusters have been variously described as "genomic islands of speciation" (Turner et al., 2005) or "genomic islands of divergence" (Ellegren et al., 2012;Jones et al., 2012;Kronforst et al., 2013;Nosil et al., 2009;Poelstra et al., 2014;Shaw & Mullen, 2011). These regions are thought to arise either by selection alone or due to inversions (Coyne & Orr, 2004;Cruickshank & Hahn, 2014;McGaugh & Noor, 2012;Moyle & Payseur, 2009;Nosil & Feder, 2012;Orr et al., 2004;Presgraves, 2010). We did not observe any evidence of inversions in our analysis and a parsimonious interpretation of our results is that selection has produced the F I G U R E 4 The upper circos plot displays fixed SNPs, G ′′ ST , Tajima's D, and linkage disequilibrium (r 2 ). Tick marks are in units of megabases. Note that the scale of the microchromosomes has been increased for clarity. Regions of increased differentiation (blocks falling within the top 1% of G ′′ ST outliers) are highlighted with gray boxes. The lower panels show increasing degrees of magnification of the divergent regions. Note that the regions containing cd36 and mlph are highlighted in yellow and blue, respectively. Tajima's D in A. m. marmoratus and A. m. speciosus  divergent regions. However, structural variants, and inversions in particular, are difficult to identify with short-read sequencing, so this conclusion will require further exploration (Lledó & Cáceres, 2013).

| Genes in divergent regions
Because the divergent regions contained clusters of fixed SNVs, we identified 5 kbp blocks that contained at least one fixed SNV. This reduced the number of genes to 97. These blocks tended to cluster together in the center of the five divergent regions (Figure 4). This set of genes contained two potential pigmentation genes: mlph and cd36. Eight SNPs were within mlph introns and five within cd36 introns. While the exons of both genes have coverage across all bases, in neither gene did a fixed SNP intersect a coding sequence. If either of these genes is contributing to the differences in pigmentation between the subspecies, this observation suggests that selection may be acting on regulation of gene expression. Additional research is needed to determine how coloration is being regulated at the gene level.
Neither mlph nor cd36 are among known genes in the pteridine pigmentation pathway. Mlph regulates melanosome transport within melanocytes (Kuroda & Fukuda, 2004;Provance Jr. et al., 2002;Schluter, 2001;Wu et al., 2002). In humans and mice, mutations in mlph are characterized by reduced pigmentation in skin and hair (Matesic et al., 2001). However, phenotypic divergence between A. m. marmoratus and A. m. speciosus does not appear to involve differences in melanin pigmentation. Cluster of differentiation 36 is a transmembrane glycoprotein that helps direct the uptake of fatty acids (van Bennekum et al., 2005). In mice and silkworms, cd36 appears to regulate the uptake of carotenoids (Doebeli & Dieckmann, 2003;Sakudoh et al., 2010Sakudoh et al., , 2013van Bennekum et al., 2005). Thus, cd36 appears to be an important candidate locus to explain differences in carotenoid pigmentation.

| HPLC
Because carotenoid pigments are known to be honest indicators of quality (Kondrashov & Kondrashov, 1999;Olson & Owens, 1998;Schluter, 2009), because cd36 was a top hit in our genome scan, and because carotenoids usually produce orange or yellow pigmentation, we used HPLC to characterize the type and quantity of carotenoid pigments in the orange-pigmented and non-orange-pigmented skin of the dorsolateral head of the two subspecies. We found that esterified-xanthophyll esters were significantly more abundant in the conspicuous orange-pigmented skin of A. m. marmoratus. While we did not measure pteridine pigments, the fact that we did not observe any genes associated with the pteridine biosynthetic pathway in our genome scan suggests that the pathway is not under selection in the two subspecies. Additional studies are necessary to determine whether cd36 is differentially expressed in the epidermal tissue of A. m. marmoratus, but the HPLC analysis strongly suggests that differences in carotenoid abundance explain the differences in phenotype.

| Sex chromosome divergence
When we investigated the divergence on macrochromosomes and microchromosomes, we observed that microchromosome LGb was significantly more divergent between the subspecies. This divergence is similar to the high divergence observed on male sex chromosomes in other pairs of incipient species such as Ficedula flycatchers (Ellegren et al., 2012). However, in A. carolinensis LGb is the female sex chromosome (Alföldi et al., 2011;Losos, 2009). Assuming that LGb also is the female sex chromosome in A. marmoratus, its divergence could have at least three explanations: (1) reduced effective population size of the female chromosome is increasing neutral divergence (Charlesworth et al., 1987;Williams, 1972); (2) selection on female genes in LGb is driving divergence (Charlesworth et al., 1987;Nicholson et al., 2007;Ord & Martins, 2006;Vanhooydonck et al., 2009); or (3) limited female dispersal which would result in increased divergence on the female chromosome.
Sex-biased dispersal has been observed in other Lesser Antillean anoles (Johansson et al., 2008;Losos, 2009). Phylogenies generated from maternally inherited mitochondrial DNA in these subspecies form monophyletic groups (Schneider, 1996), which suggests that female dispersal may be limited. Population genetic theory suggests that male sex chromosomes should diverge even more rapidly because they do not recombine and are subject to more cell divisions (Graves, 2006;Maan & Seehausen, 2011). Thus, an additional possibility is that LGb may actually be the male sex chromosome in A. m. marmoratus.

| CON CLUS ION
The discovery that a gene (cd36) potentially involved in carotenoid sequestration is under divergent selection in A. m. marmoratus and A.
m. speciosus, and that the conspicuously pigmented skin along the dorsolateral head of adult male A. m. marmoratus contains significantly more carotenoid pigments, suggests that cd36 may be contributing to the differences in pigmentation between these two subspecies.
This is particularly exciting because carotenoid pigments play important roles in sexually selected social signals and are considered to be honest indicators of quality (Olson & Owens, 1998;Olsson et al., 2013). Despite research on the genetic basis of carotenoidbased phenotypes in reptiles and fish (García-de Blas et al., 2013;Goodwin, 1984;Macedonia et al., 2000;Olson & Owens, 1998;Pike et al., 2010;Tripathi et al., 2009;Walsh et al., 2012;Ziegler et al., 2000), no carotenoid genes have been previously identified to be the targets of natural selection. Because cd36 regulates the deposition of differently colored carotenoid pigments in mice and insects, our results suggest that it may be a key gene regulating the deposition of carotenoid pigments in animals. Gene expression and immunohistological analysis of cd36 and mlph will help further characterize whether these genes are active in the pigmented epidermis of A. m. marmoratus and A. m. speciosus.
From a broader perspective, sexual selection and the incidence of sexual dichromatism have been shown to correlate with increased species diversity (Barraclough et al., 1995;Saks et al., 2003). Anolis lizards are remarkably diverse, with more than 400 described species and sympatric species almost always have different colored, conspicuous dewlaps (Losos, 2009;Nicholson et al., 2007;Williams & Rand, 1977). Within anoles, differences in coloration appear to evolve prior to the evolution of larger morphological changes, suggesting that divergence in pigmentation is an important component of speciation in anoles (Endler, 1980(Endler, , 1983Losos, 2009). The regions of clustered divergence observed in our study and the genes contained within them provide insight into the early stages of speciation within anoles. Future work in similar anoline systems such as the A. distichus, brevirostris, and apletophallus species complexes will help characterize whether the same regions and genes are contributing to early stages of speciation in anoles (Lambert et al., 2013;Ng & Glor, 2011;Olson & Owens, 1998;Stapley et al., 2011). Additional

ACK N OWLED G M ENTS
We thank Martha Muñoz for help organizing the collection of many of the specimens as well as our enthusiastic field assistants Juanita